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Collisionless suspensions of inertial particles (finite-size impurities) are studied in two- 
and three-dimensional spatially smooth flows. Tools borrowed from the study of random 
dynamical systems are used to identify and to characterise in full generality the mecha- 
nisms leading to the formation of strong inhomogeneities in the particle concentration. 

Phenomenological arguments are used to show that in two dimensions, the posi- 
tions of heavy particles form dynamical fractal clusters when their Stokes number (non- 
dimensional viscous friction time) is below some critical value. Numerical simulations 
provide strong evidence for the presence of this threshold in both two and three di- 
mensions and for particles not only heavier but also lighter than the carrier fluid. In two 
dimensions, light particles are found to cluster at discrete (time-dependent) positions and 
velocities in some range of the dynamical parameters (the Stokes number and the mass 
density ratio between fluid and particles). This regime is absent in the three-dimensional 
case for which evidence is that the (Hausdorff) fractal dimension of clusters in phase 
space (position-velocity) remains always above two. 

After relaxation of transients, the phase-space density of particles becomes a singu- 
lar random measure with non-trivial multiscaling properties, whose exponents cannot be 
predicted by dimensional analysis. Theoretical results about the projection of fractal sets 
are used to relate the distribution in phase space to the distribution of the particle posi- 
tions. Multifractality in phase space implies also multiscaling of the spatial distribution 
of the mass of particles. Two-dimensional simulations, using simple random flows and 
heavy particles, allow the accurate determination of the scaling exponents: anomalous 
deviations from self-similar scaling are already observed for Stokes numbers as small as 
10-*. 



1. Introduction 

Suspensions of dust, droplets, bubbles or other kind of small particles in turbulent 
incompressible flows are present in many stirring and mixing problems encountered in 
both natural and industrial situations. Such impurities typically have a finite size and a 
mass density different from that of the carrier fluid. They cannot be described as simple 
passive tracers, that is point-like particles with negligible mass advected by the fluid; an 
accurate model for their motion must take into account inertia effects due to the finiteness 
of their sizes and masses. These inertial particles interact with the fluid through a viscous 
Stokes drag and thus their motion typically lags behind that of passive tracer particles. 
The dynamics of the latter is governed by a conservative dynamical system when the 
carrier flow is incompressible (because volume is conserved), but inertial particles have, 
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as we shall see, dissipative dynamics. While an initially uniform distribution of passive 
tracers remains uniform at any later time, the spatial distribution of inertial particles 
develops strong inhomogeneities. 

Such a phenomenon, frequently referred to as preferential coneentration, is associ- 
ated with the presence of regions with either extremely high or low concentrations. 
Their characterisation plays an essential role in the understanding of natural and in- 
dustrial phenome na. Instances are op timisation of combustion processes in the design 



of Diesel engines dElperin et a?Jl2003l). the g rowth of rain drops in subtropical clouds 
ijPinskv fc Khainlll997t iFalkovich et a/]l2002|) . the a ppearance of planetes i mals in the 
early stage of planet formation in the Solar system | Weidenschilling' ' 1 9 9 5'; 'Cuzzi et al\ 
2^^ ) , coexistence proble ms between several species of plankton (.Sauii'es &: Yamazakil 
1995 : iKkrolvi et oi.ll200fll) and other environmental problems (see, e.g . Iseinfeldlll98fil) . 



For such applications it is recognised that a key problem is the prediction of the collision 
or reaction rates and t heir associated typical time scales. The traditional way to estimate 
the latter goes back to lTavlod l)l92lj) and makes use of diffusion theory. The time scales 
obtained in such approaches generally exceed by one or several orders of magnitude those 
observed in experiments or numerical simulations (see, e.g., [Su ndar a m fc C ollins 1997). 
A full understanding of particle clustering and, in particular, of the fine structures ap- 
pearing in the mass distribution is crucial for identifying and quantifying the mechanisms 
responsible for the drastic reduction in time scales. 

We propose in this paper an original approach leading to a systematic description of 
inertial particles clustering. This approach is in part inspired by recent major break- 
throughs in the study of passive s calar advection by tu rbulent flows, using Lagrangian 
techniques (see, e.g., the review of IFalkovich et oiJl200lh . Preferential concentrations are 
due to the convergence of inertial particle trajectories onto certain sets in the position- 
velocity phase space called attractors, which are usually fractals. Of course, these sets 
are dynamically evolving objects and depend on the history of the carrier flow. When 
interested only in the spatial distribution of particles, one clearly has to consider the 
projection on positions of the phase-space attractor. Projection of singular sets are them- 
selves generally singular. As a consequence, an initially uniform distribution of particles 
will tend to become singular at large times, after relaxation of transients. This is the 
basic mechanism responsible for the formation of clusters of particles. Use of dissipative 
dynamical systems tools and, in particular, of methods borrowed from the study of ran- 
dom dynamics, allows a rather complete characterisation of the particle distribution at 
those scales where the carrier flow can be considered smooth in space. As we shall see, 
the statistical properties of the spatial distribution of particles in these clusters can be 
characterised as a function of the dynamical parameters (which depend on the physical 
properties of the fluid and on the masses and sizes of the particles) . 

We focus on a simple model which captures most qualitative aspects of inertial parti- 
cles. In this simplified dynamics, only two effects are included: the Stokes drag propor- 
tional to the velocity difference between the particle and the fluid and an added mass 
term proportional to the acceleration of the fluid volume displaced by the particle. These 
two terms are associated to two parameters: the Stokes number St proportional to the 
cross section of the particle and the mass density ratio between the particle and the 
surrounding fluid. We identify and characterise the different regimes of the dynamics 
appearing when varying these two parameters. 

The paper is organised as follows. In |21'W6 derive a simple model for the dynamics 
as an approximation of the full Newton equation describing the motion of particles; in 
particular, we justify its relevance in some asymptotic range of the physical parameters 
of the problem. In ^ this model is used to relate the local dynamics of the particles 
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to the local properties of the carrier flow; this permits to perform a complete stability 
analysis of the particle dynamics. In fQ] after a brief introduction to the tools used in 
studies of dissipative dynamics with particular emphasis on Lyapunov exponents f H4.1|l . 
phenomenological arguments are used to extend the local approach to the global dynamics 
in order to show that below a critical value of the Stokes number the particles concentrate 
on fractal dynamical clusters in the position space f ^4.2|) : this is confirmed by numerical 
experiments presented in ti4.3l A more precise description of the statistical properties of 
these clusters is performed in fj^lwhere it is shown that the distribution of particle actually 
has multifractal properties, even at very small Stokes numbers. Finally, ^encompasses 
concluding remarks and open questions. 



2. Dynamical model for the particles 

We consider very diluted suspensions of inertial particles where collisions, particle-to- 
particle hydrodynamical interactions and retroaction of the particles on the motion of the 
fluid can be ignored. It is assumed that the carrier flow has moderate Reynolds numbers 
and that the particle radius is smaller than their dissipation scale. Under the additional 
assumption that the Reynolds number base d on the particle size a nd its relative velocity 
with respect to the fluid is sufficiently small. lMaxev fc RilevI l)l983l) proposed approaching 
the flow surrounding the small sphere by a Stokes flow. Neglecting both the effects of 
gravity and the Faxen corrections, the motion X(t) of an inertial particle is then solution 
to the following Newton equation 



WpX — nif (X, t) — 67ra/i 



X-m(X, t) 
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dt 



iu{X,t)) 



ds 



*^ ^X-m(X,s) 



ds 



(2.1) 



The dots denote time derivatives, u is the velocity field of the carrier fluid, Du/Di is its 
derivative along the path of a fluid element, is the mass of the particle, my is the 
mass of fluid displaced by the particle, a is the radius of the particle and, finally, /i and 
u are respectively the dynamic and the kinematic viscosities of the fluid. The different 
forces exerted on the particle are, in order of appearance on the right-hand side of (|2.1|) . 
the force exerted by the undisturbed flow, the Stokes viscous drag, the added mass effect 
and the Basset-Boussinesq history force. 

For rescaling purposes it is more convenient to write the equation satisfied by the 
velocity difference between the particle and the carrier fluid w{t) = X — tt(X, t). Indeed, 
the dynamics involves two velocity scales which are determined by u and w and there 
is a priori no reason why they should be of the same order of magnitude. Rescaling of 
space, time, velocity of the carrier fluid and velocity difference by the factors L, L/U, U 
and W, respectively, leads to 



dw 
'df 



13-1 d 
a dt 



(M(X,t))-/3(u;- V)m(X,<) 



5*^ 



3/3 

TT St ,/o 



ds dw 

\/t — s ds 



(2.2) 



The dynamics depends here on three parameters : the velocity ratio a = W/U, the added- 
mass factor /? = 3mf/{mf + 2mp) and the Stokes number associated to the particle 
St = {a'^U)/{3PvL). The latter can be also written as St = Re {a/Lf/{2,(3) where Re = 
UL/v is the Reynolds number of the carrier flow. The particle radius a is taken to be 
smaller than the typical (integral) scale L of the carrier flow. As a consequence, if the 
Reynolds number is not very large and the particles are not very heavy we must assume 
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that the Stokes number is much smaUer than unity. Moreover, the requirement of small 
local Reynolds number in the neighbourhood of the particle - required in Maxey & Riley 
derivation of H2.1|) - reads here Wa /v <^\.li restricts the range of admissible parameters 
to satisfying a ^ (/3 St)~^ . When both j3 and St are taken order unity, the velocity ratio 
a must be small compared to unity, meaning that the inertial-particle dynamics must be 
very close to that of ordinary fluid particles. This can only be achieved by requiring low 
values of the Stokes number. Hence, in the Maxey & Riley approach, (3 and St cannot 
be simultaneously order unity. 

We nevertheless consider here particles whose motion is described by a simplified dy- 
namics where both the second and the last term of H2.2|l are neglected. This is clearly the 
case when the typical velocity difference a, and thus the Stokes number are both very 
small. Such simplified dynamics is also relevant when the particles are much heavier than 
the fluid {(3 <^ 1), a case of importance for applications involving suspensions of liquid 
droplets or dust particles in a gas. Under one of these two hypotheses and by introducing 
what we call here the eovelocity V = X — /3 m(X, t), it is easy to see that the equation of 
motion reduces to the (2 x (i)-dimensional system 

±^ (3u(X,t) + V, (2.3a) 

V=l[(l-/3)«(X,t)-V], (2.36) 

combined with the initial conditions X(0) = Xq and V(0) = (1 — f3)u{xo,0) + aw^. 
Here Wq denotes the non-dimensional initial velocity difference between the particle and 
the flow. The dynamics described by (|2.3|) are clearly dissipative. Indeed, the divergence 
of the right-hand side of H2.3I) with respect to the phase-space variables {x,v) reduces 
in the case of divergence-free carrier flow to —d/St, which is negative. This implies 
that all volumes in position-covelocity phase space are uniformly contracted during time 
evolution. 



3. Local analysis of the dynamics 

As a first step in describing the mechanisms responsible of particle clustering, we 
study the time evolution of the position-covelocity phase-space separation R = (5X, 5V) 
between two infinitesimally close trajectories. Since the sizes of inertial particles are 
within the dissipation range of the carrier flow, the velocity field u can be considered 
spatially smooth. Velocity gradients are then bounded and velocity differences between 
two neighbouring locations are linearly proportional to their separation. The separation 
R(f) thus obeys an equation obtained by linearising H2.3|l . namely 



R = XtR, 



Mt = 



Pa{t) 

1-/3 



St 



a{t) 



1 , 

Yt' 



(3.1) 



Here <t denotes the strain matrix of the carrier flow along the trajectory of a particle: 
(Jij{t) = djUi{X{t),t), and Xd stands for the d-dimensional identity matrix. Qualitative 
insight into the local dynamics of the particles is provided by making the stability analysis 
of H3.1(l . The eigenvalues of the evolution matrix Mt depend on the local structure of the 
carrier flow, as we now show. It is easily checked that A^t is a solution to the second-order 
equation 



St 



(3.2) 
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Figure 1 . Local structure of the carrier flow in three dimensions ; four cases are distinguished 
accordingly to the nature of the three eigenvalues of the strain matrix: if all three are real 
(hyperbolic case), either two of them are negative and one is positive (o) or two of them are 
positive and the third is negative (6). When two eigenvalues are complex conjugate (elliptic 
case), the third one is real and can be either negative (c) or positive (d). 



where is the 2d x 2d block-diagonal matrix given by 

(3.3) 



(Tit) Od 
Od T{t) 



and Od is the d-dimensional square matrix with all entries zero. A straightforward conse- 
quence is that, to a given eigenvalue 7 of the strain matrix, correspond two eigenvalues 
of the evolution matrix Mt solutions of the quadratic equation 

.^+(]^-/37)x-|-0. (3.4) 

The stability analysis of the evolution matrix requires to express the sign of the real part 
of the solutions to this second-order equation as a function of 7, and thus depends on 
the local structure of the carrier flow. 

Let us examine the conse quence, first in two dimens ions. Depending on the sign of the 
Okubo- Weiss parameter Q (*Okubol ll97(it IWeis jHool defined as the determinant of the 
strain matrix cr and thus satisfying 7^ = Q, the incompressible flow of the carrier fluid 
separates the space in regions of two kinds: the hyperbolic regions where Q > and where 
the strain matrix cr has two opposite real eigenvalues and the elliptic regions where Q < 0, 
so that the two eigenvalues are imaginary numbers. A simple analysis of (|3.4() shows that 
in the hyperbolic regions, the evolution matrix Mt of the particle dynamics has three 
stable and one unstable eigendirections in the four-dimensional position-covelocity phase 
space. In the elliptic regions, one must distinguish two cases: if /3 < 1 (particles heavier 
than the fluid), there are two stable and two unstable eigendirections, whereas if /3 > 1 
(lighter particles), there are always four stable directions. This local analysis suggests 
that in two dimensions, heavy particles are excluded from the vortices (elliptic regions) 
and tend to concentrate within the fllaments (hyperbolic regions) of the carrier flow, 
whereas light particles tend to cluster in the vortices. 

In three dimensions, there are four possible local structures of the incompressible carrier 
flow represented schematically in figure They can also be divided into elliptic and 
hyperbolic cases depending on the presence or not of complex eigenvalues. 

• Hyperbolic regions (a, b): the local dynamics of the particles has five stable and one 
unstable eigendirections in the case (a) and four stable and two unstable in the case (&); 

• Elliptic regions (c, d): there is generally at least one unstable eigendirection; the 
exception is case (c) for light particles {f3 > 1); all six eigendirections are then stable 
provided the real negative eigenvalue 7 of the strain matrix satisfies 7 < —St/ (3. 

The three-dimensional situation looks hence similar to the two-dimensional one: the 
local analysis suggests that particles lighter than the carrier fluid concentrate in the 
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rotation regions while heavier particles tend to be ejected from them and concentrate 
in the strain regions. The model is hence able to catch a qualitative property of inertial 
particles dynamics which plays a central role in the pheno menological understandi ng of 



their preferential concentration. As stated for instance in ISauires fc EatonI ||1991|) . the 



appearance of inhomogeneities in the distribution of particles is generally explained by 
the presence of persistent structures in the carrier flow. The latter are responsible for 
a deterministic motion of the particles during which their inertia affects their relative 
motion with respect to the carrier flow: phenomenology indicates that heavy particles are 
ejected from eddies while light particles tend to concentrate in their cores. This is indeed 
what our local analysis predicts. However we shall see in the next section that, although 
preferential concentrations of inertial particles are triggered by the local mechanisms just 
discussed, it is the dissipative character of the dynamics which is solely responsible for 
the eventual strong clustering of particles. 

4. Threshold in Stokes number for the formation of fractal clusters 

4.1. Dissipative dynamics and Lyapunov exponents 

The temporal evolution of the separation R(t) between two infinitesimally close trajec- 
tories of the phase space can be expressed in terms of the Green function J7t associated 
to the linearised dynamics H^H|I . We indeed have 



where Texp denotes the time-ordered exponential of matrices and A4t is the evolu- 
tion matrix defined in 13. 1() . The symmetrical matrix Jt Jt has positive eigenvalues 
that can be written e^^"^*-*^ e^^^''*-*-' *. The [ijS are called the stretehing rates (or 
the finite-time Lyapunov exponents) and are generally labelled in non-increasing order 
A*i(i) ^ • ■ • 5^ M2ti(^)- They measure the time-evolution of inflnitesimal volumes of phase 
space. ^1 measures the exponential growth rate of the distance between two neighbour- 
ing trajectories, fXi + fi2 measures that of areas defined by three trajectories, etc. The 
sum of the 2d stretching rates controls the time evolution of phase-space 2(i-dimensional 
volumes and is expressible in terms of the trace of the Green matrix Jt . It is easily shown 
that + ■ • ■ + /i2d = —d/ St. 

The long-time behaviour of the local dynamics is dominated by the almost-sure con- 
vergence of the stretching rates to t he classical Lya punov exponents Xj = limt^oo (t) ■ 
The multiplicative ergodic theorem l|Oseledetslll968|) ensures that, under some ergodicity 
hypothesis on the dynamics, the Lyapunov exponents are independent of both the reali- 
sation of the random carrier flow and of the peculiar trajectory X(t) around which the 
linearised dynamics is considered. The Lyapunov exponents are linked to many funda- 
mental features of the dynamics. For instance, when the largest Lyapunov exponent Ai 
is negative, the stability of the linearised system is ensured and all the particle trajecto- 
ries are converging together, thereby leading to a somewhat degenerate statistical steady 
state in which all the mass is concentrated in discrete time-dependent points in phase 
space. However when Ai is positive, the situation is more complex: the dynamics is then 
said to be chaotic and stability is ensured only if the initial phase-space separation R(0) 
is orthogonal to the subspace generated by the eigendircctions associated to the positive 
Lyapunov exponents. 

When Ai > the long-time dynamics has richer features, characterised by the conver- 
gence of the particle trajectories to complex dynamical structures called attractors. These 







(4.1) 




Figure 2. The Lyapunov dimension is defined by interpolating to non-integer values the 
exponential growth rate of phase-space fc- dimensional volumes obtained from the sum of the k 
largest Lyapunov exponents. 



are generally (evolving) fractal sets of the phase space and can be characterised by their 
HausdorfF dimension dn^ which is expected to coincide with the box-counting dimension 
in non-degenerate cases. The appearance of attractors can be seen as the result of a 
competition between stretching and folding effects that occur during the chaotic motion 
of the particles. As a consequence, the properties of these attractors, and in particular 
their dimension, depend on the stretching rates of the dynamics. The positive //j 's are 
responsible for stretching in their associated eigendirections, while the negative rates give 
contra ction and hence lead to folding. Following this kind of approach, Kaplan fc Yorkj 
l)l979() proposed to estimate the dimension of the attractor by the Lyapunov dimension, 
defined as 

= J- ^i-t— ^i— where Ai H h A,/ ^ and Aj H + Aj+i < 0. (4.2) 

A, 7+1 

This non-random number can be interpreted heuristically as the dimension of phase- 
space objects that keep a constant volume during time evolution. As already stated 
above, the long-time evolution of fc-dimensional volumes is governed by the sum of the 
k first Lyapunov exponents. The Lyapunov dimension is obtained by interpolating this 
su m to non-integer values and determining where it vanishes (see figure It was shown 
bv lDouadv fc Pester id l)l980(l that it serves as a rigorous upper-bound for the Hausdorff 
dimension djy; this will be used in next subsection to show that below a critical value 
of the Stokes number, the particles form fractal clusters in the position space. As we 
shall see in |S1 the Lyapunov dimension is actually equal to the information dimension 
associated to the steady-state phase-space density of particles and is thus related to the 
small-scale properties of the mass distribution of inertial particles. 

4.2. A heuristic approach in two dimensions 

We now give a heuristic argument which predicts the presence of a threshold in St okes 
numb er for clustering of heavy particles. This result, already described briefly in iBed 
l)2003j) . is here presented in more details. 

As already stated, the particles concentrate on dynamical attractors in the position- 
covelocity phase space. The main idea of the argument presented here comes from a trivial 
observation. Focusing on positions while ignoring covelocities amounts to projecting the 
attractor from the 2d-dimensional phase space onto the rf-dimensiona l space of parti cle 
positions. A standard result on the geometry of fractal sets (see, e.g.. lFalconeitil986|) is 
that, when the Hausdorff dimension dn of the attractor is less than the dimension of 
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the projection space (here, d), the projection of the fractal is itself a fractal set with 
Hausdorff dimension dn- If however dn > d, the projection has dimension d. Now, using 
the Lyapunov dimension d^ as a rigorous upper bound of the fractal dimension dn of 
the phase-space attractor, the condition < d is a sufficient condition for particles to 
form fractal clusters in position space. The next step is to find a condition ensuring that 
the Lyapunov dimension is less than the space dimension. From its definition (|4.2I) . it is 
clear that d^, < d if the sum of the d largest Lyapunov exponents is negative. Hence, the 
problem reduces to finding a sufhcient condition for Ai + • • • + < 0. 

To estimate the d largest Lyapunov exponents we first consider the stability exponents 
defined as 

= lim - hiaj(t), (4-3) 

t — *oo t 

where the aj (t) 's are the eigenvalues of the Jacobi matrix J't defined in (|4.1|) , which are 
usually complex numbers. Although there is no equivalent of the Oseledets ergodic the- 
orem for such complex exponents, it is usually supposed that the limit in (|4.3() exists in 
generic situations and depends neither on the realisation of the carrier velocity field, nor 
on th e particular trajectory along which they are calculated (see, e.g. Goldhir sch et all 
119871^ . The real parts of the stability exponents generally differ from the Lyapunov expo- 
nents. It is however possible to obtain inequalities between these two sets of exponents, 
using Browne's theorem (see, e.g., for details M arcus fc Mine 1992). We will also need 
the strain exponents r]i analogous to the /ij's but this time associated to the {d x d)- 
dimensional matrix Texp dr (t{t). 

We now focus on the case of heavy particles {(3 < 1) embedded in a two-dimensional 
incompressible carrier flow. Note first that, as we have seen from the local analysis in the 
previous section, the heavy particles tend to concentrate in the filamentous structures 
of the flow (hyperbolic regions). The eigenvalues associated to the local dynamics of the 
particles are given by the quadratic relation H3.4I) in terms of those associated to the 
strain matrix of the carrier flow. Now comes the heuristic point. Let us assume that the 
same quadratic relation allows us to express, at least in an approximate way, the stability 
exponents a^'s as a function of the strain exponents ry^'s defined abovef. The latter are 
expected to be real since they are dominated by the hyperbolic regions of the flow. In 
two dimensions there are only two strain exponents +1] and —rj. Under these hypotheses 
a sufficient condition for fractal clustering of the particles can be written 

St^^[p-2 + 2VI-/3 + /32) . (4.4) 

Observe that the (non-dimensional) strain exponent rj is calculated along the path of 
inertial particles and not along that of fluid particles. As the heavy particles spend a 
long time in the hyperbolic regions, the strain should satisfy rj ^ \fL/U, where A/ 
denotes the positive Lyapunov exponents associated to the dynamics of fluid particles. 
Using this relation in (|4.4I) we obtain that there is clustering if 

^^"^W^ (/3-2 + 2Vl-/? + /32). (4.5) 

Note that this sufficient condition gives heuristic evidence that clustering already occurs 
at small Stokes numbers, although the Stokes drag is then very strong, and thus the 
inertial particle motion is very close to that of simple passive tracers. This is consistent 

f This assumption would be exact if one could replace time-ordered exponentials appearing 
in the definition of the various Green functions by ordinary exponentials; of course, this is not 
the case here. 
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with predic tions ma de _by IBalkovskv et a.l\ l|200ll) . Their approach is based on the ob- 
servation of iMaxevI l|l987l) that the discrepancy between inertial particle dynamics and 
passive tracers occuring when St <^ \ can be captured by a spatial Taylor expansion 
of the particle velocity field, leading to a model in which the particle is advected by a 
synthetic flow comprising a small compressible component. Such an expansion allowed 
Balkovsky et al. to show that moments of particle density grow exponentially in times, 
which is indeed indicative of the formation of clusters. 

The existence of a threshold in Stokes number for particle clustering requires that 
for large-enough values there be no clustering. It is easily checked from equations H2.3|l 
governing inertial-particle dynamics that, when increasing the Stokes number, the particle 
are less and less influenced by the carrier flow. More precisely, in the limit St ^ oo, their 
motion becomes purely ballistic and their distribution in position space tends to become 
uniform, ensuring absence of fractal clusters for large-enough Stokes number. 

Note that the same type of arguments when applied to light particles (with (3 > 1) 
whose dynamics is dominated by elliptic regions of the flow, would predict that they 
always concentrate on points located in the cores of vortices (clusters of dimension dn = 
0). As we shall see from the numerical experiments of the next subsection, there are 
ranges of values of the parameters (3 and St for which this strong clustering occurs in two 
dimensions, but for other values the light particles can also form either fractal clusters or 
fill the whole space. This indicates the limits of our heuristic approach. The arguments 
used in the case of heavy particles do not apply to light particles in the elliptic regions 
since the eigenvalues of the local dynamics are then complex and cannot be used to 
estimate the stability exponents. Moreover, in order to match the limit St in which 
the motion of light particles recovers that of fluid particles, it is clear that their dynamics 
is influenced by their stay in the hyperbolic regions, even if this stay is shortened by 
inertial effects. The phenomenological derivation of the sufficient condition H4.5|l was just 
meant to give an idea of the physical mechanisms responsible for particles clustering. Its 
qualitative validity is however confirmed numerically, as we shall now see. 

4.3. Numerical evidence for fractal clustering in two and three dimensions 

In our simulations the inertial particles positions are confined to a periodic domain of 
size L. The prescribed carrier flow is a solution to the forced Stokes equation 

dtU = i^Au + f{x,t), V-M = 0, (4.6) 

where f{x,t) is a space-periodic, isotropic and homogeneous Gaussian random forcing 
concentrated at spatial scales of the order of the box size L and delta-correlated in time. 
Such a carrier flow provides control over the statistical properties of the Gaussian velocity 
field u and allows accurate numerics. In practice we take v = 1 and assume that the 
velocity field has the lowest reasonable number of spatial Fourier modes for which isotropy 
is approximately ensured: 9 in two dimensions and 27 in three dimensions! ■ Note that 
two-dimensional flows satisfying (|4.6(l were also considered by ISig^irgeirsson fc StuartI 
who proved in this setting existence of the random dynamical attractor for very 
heavy particles. 

The individual trajectories of inertial particles are integrated using a fourth-order 
Runge-Kutta scheme. Choosing a carrier fluid velocity with rather few Fourier modes 
makes it unnecessary to interpolate the velocity at particle locations, since we can calcu- 
late it by direct summation of the Fourier series. Fully resolving the finest structures of 



f Fewer Fourier modes may actually lead to degenerate situations and in particular to insuf- 
ficiently mixing dynamics. 
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(a)d^2 (b)d^3 

Figure 3. Phase diagrams in the parameter space St) for d = 2 (a) and d — 3 (b) represent- 
ing the different regimes of the dynamics corresponding to different behaviours of the Lyapunov 
exponents. For parameter values in the white area (such that Ai + • ■ ■ + Ad > 0), the distribution 
of particle fills the whole position space. In the light-gray area the Lyapunov dimension is less 
than d (i.e. Ai -I- ■•■-)- Ad < 0) and particles concentrate on fractal clusters. The darker gray 
area which is only present in the two-dimensional case corresponds to the regime for which all 
Lyapunov exponents are negative: the asymptotic distribution of particles is then atomic. 



the spatial particle distribution requires a strict control over the small-scale dynamics, 
as is here the case. 

We present in this section measurements of the Lyapunov exponents associated to the 
dynamics of the particles. To est imate these quantitie s, we use the sta,ndard method de- 
veloped bv lBenettin et all l|l98(t) - see also the book bv lCrisanti et all ()l993tl for a precise 
description of this method. The idea is to integrate the time evolution of 2d infinitesimal 
separations i?i , . . . , R2d governed by the linearised dynamics H3.1(l and along the path of 
a given particle. The Lyapunov exponents are obtained from the exponential growth rates 
of the distances, surfaces, volumes, etc. defined by these infinitesimal vectors. In order 
to prevent numerical errors from accumulating and the length of the Rj 's from increas- 
ing too rapidly, their time evolution must be accompanied by frequent renormalisation,f 
using, for instance, the Gram-Schmidt procedure. In practice, to obtain the Lyapunov 
exponents for given values of the parameters f3 and St, we integrate the linearised system 
for 10^ turn-over times performing a renormalisation of the vectors Ri , . . . , i?2d every 
tenth of turn-over time. 

Figure 13 (a) shows the phase diagram obtained numerically in the two-dimensional 
case, which divides the parameter space (/3, St) into three different regimes. The light- 
gray area represents those values for which the sum of the two first Lyapunov exponents 
is negative. The Lyapunov dimension is there smaller than the dimension d = 2 of the 
position space and the particles form fractal clusters. The second regime appearing in 
the dynamics is when the sum of the two first Lyapunov exponents is positive (white 
area) . This corresponds to the case when the Lyapunov dimension is larger than two and 
the particles fill the whole domain. The last regime, represented by the dark-gray area 
in figure 121 (a), corresponds to values of the dynamical parameters for which the largest 
Lyapunov exponent Ai is negative. This strong clustering happens only when the particles 
are lighter than the fiuid. They form there point-like clusters and the associated mass 

f Here the term "renormalisation" is to be taken literally, meaning just "making the norm 
finite" . 
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distribution is said to be atomic. In the case of bounded domains (such as considered 
here) it corresponds to the large-time convergence of ah particle trajectories towards a 
single trajectory. 

In three dimensions, the situation is rather different. As shown in figure|31 (6), there 
are only two different possible regimes in the particle dynamics: one corresponding to 
a Lyapunov dimension larger then the position-space dimension and in which particles 
fill the whole space; the other one associated to a Lyapunov dimension smaller than d 
and where particles concentrate onto fractal clusters. Surprisingly, there is no evidence 
for the presence of a region in the parameter plane where the particles cluster on a set 
of dimension less than two. A systematic investigation for a large number of different 
values of the parameters /3 and St was done. In particular, the Lyapunov exponents were 
calculated for forty different values of St aX (3 = 2.9 for which clustering would have been 
expected. Results indicate that the sum of the two largest Lyapunov exponents always 
remains positive. This observation contradicts the phenomenological understanding of 
clustering which predicts that light particles may form clusters of dimension between 
one and two in the core of the vortices. The absence of such a behaviour is not due to the 
choice of the model, otherwise point-like clusters would also be absent in two dimensions. 
A possible explanation is related to the properties of the carrier flow considered here. It 
could be that in our setting, the velocity field u does not present rotational structures 
which are strong enough and sufficiently persistent to permit such a clustering of the 
particles. On the one hand, we have seen from the local analysis of f|31 that the elliptic 
regions of the carrier flow are stable for light particles only if the Stokes number is 
sufficiently small. On the other hand, St has to be sufficiently large for the particles to 
have enough inertia to separate from the fluid motion. It may happen that these two 
conditions are here incompatible. It would be of interest to confirm the absence of point- 
like clusters by systematic investigations in different flow fields (and in particular realistic 
ones); this will be done in future work. 

The separatrix between the regimes associated to different signs of the sum of the d 
largest Lyapunov exponents has a discontinuity at /3 = 1 in both two and three dimen- 
sions. This singular behaviour is due to the degenerate dynamics associated to particles 
with the same density as the fluid (such particles are usually referred to as neutrally 
buoyant particles). It is indeed clear from H2.3|l that when /3 = 1, the covelocity does not 
depend anymore on the velocity of the carrier flow; it simply relaxes exponentially from 
its initial value to zero. As a consequence, the subspace of covelocities is an eigenspace 
for the linearised dynamics associated to a Lyapunov exponent equal to —1/ St and with 
multiplicity d. Incompressibility of the carrier flow then implies that Ai + • • • + = on 
the whole line /3 = 1 of the parameter space, which partly explains the singular behaviour 
observed here. 

More quantitative information is provided by the behaviour of the Lyapunov dimension 
as a function of the Stokes number. Figure 01 shows the results obtained from both 
two- and three-dimensional simulations in the case of very heavy particles {(3 — 0). First, 
we observe that the dimension of the attractor indeed tends to the space dimension d at 
low Stokes numbers. More precisely, in both two and three dimensions, the Lyapunov di- 
mension behaves in this asymptotic regime as c^l — d—C St^ . Such a quadratic beh aviou r 
near vanishing Stokes numbers has been obtained theoretically bv lBalkovskv et ali lj200lh 
using the method of advection by a synthetic compressible flow, already mentioned. At 
large Stokes numbers, the Lyapunov dimension becomes larger and larger, eventually 
reaching 2d when St ^ oo. Between these two asymptotic regimes, there is a range of 
Stokes numbers for which the Lyapunov dimension - and thus the fractal dimension of the 
attractor - is smaller than d. For such Stokes numbers, the particles form fractal clusters. 
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Figure 4. Difference cLl — d between the Lyapunov dimension associated to the dynamics 
of very heavy particles (/3 = 0) and the position-space dimension, as a function of the Stokes 
number St (circles: d = 2, squares: d = 3). The critical Stokes number corresponds to the value 
for which cLl = d. Upper-left inset : same in log-log coordinates showing evidence for a quadratic 
behaviour at low Stokes numbers. 

The graphs of the Lyapunov dimension against the Stokes number, shown in figure 01 
display a minimum around St — 0.09 for c? = 2 and St — 0.14 for d — i. These are the 
Stokes numbers for which the strongest clustering takes place. Finally, these experiments 
also confirm the existence of a critical value for the Stokes number below which fractal 
clustering occurs; on the figure it corresponds to the value for which the curves cross the 
Hne d = dL (around St = 0.24 for d = 2 and St = 0.32 for d = 3). 

Qualitative insight into clustering is provided by figure where snapshots of the po- 
sition of 10^ very heavy particles illustrate the different regimes arising in the two- 
dimensional dynamics when increasing the Stokes number. First, when the Stokes number 
is very small (here St = 10"'^ for figure a), the particles are rather densely distributed 
but it is possible to observe fractal clusters with a dimension very close to the space 
dimension d. This picture illustrates the highly singular nature of the vanishing Stokes 
number limit. Note also that, as the Stokes number becomes very small, the dynamical 
fractal clusters are becoming denser and more space-filling. Conversely, as we increase St, 
the dimension of the clusters decreases and large empty areas appear (figure |31 h) . This 
process continues up to the value of St corresponding to the strongest clustering. Above 
this value, the dimension of clusters increases as a function of St and caustic-like struc- 
tures become conspicuous in the particle spatial distribution (figure Elc). Finally, when 
the Stokes number is larger than the critical value, the particles fill the whole domain, 
but in a nonuniform way (figure |31 d). 

5. Statistical properties of the mass distribution 

5.1. Phase- spaee density fluctuations and anomalous sealing 

In the previous section the Lyapunov dimension was used as an estimate to determine 
the presence or not of fractal clusters of particles. This dimension is not only an upper 
bound on the fractal dimension of the dynamical attractor; it also provides additional 
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Figure 5. Snapshots of the positions oi N = 10^ heavy particles (/3 — 0) associated to four 
different Stokes numbers as labeUed. (a), (b) and (c) correspond to values smaller than the 
threshold, so that particles form fractal clusters, (d) corresponds to a Stokes number larger 
than the critical value, so that the particles fill the whole domain. 

information on the small-scale properties of the particle distribution. Consider an initial 
density of particles in the phase space which is uniform when projected on the position 
space and may or may not possess some dispersion in covelocities at each position. At 
a given time t sufficiently large to be close to the asymptotic regime, the phase-space 
density field f{x, v, t) is generally singular: its support is exactly the random dynamical 
attractor toward which all individual trajectories converge. This limiting density is a 
random equivalent of the Sinai-Ruelle-Bo w'en (SRB) measure, well-known in the study 
of dissipative dynamical systems (see, e.g.. lEckmann fc RuelldllQSSt lYoundl2002|) . The 
density distribution along the attractor can be characterised by the mass of particles 
contained in small balls of radius r centred on the attractor. The requirement to be on 
the attractor is automatically satisfied if we let the ball centre follow a given trajectory 
(X(t), V(t)) in phase space. Let us thus introduce the phase-space mass in a phase-space 
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ball or radius r around the trajectory: 



TO^(X(t),V(t),i) = / dxAv f{y^{t)+x,Y{t)+v,t). (5.1) 

Note that x and v are here non-dimensional variables, so that it is meaningful to use 
the Euclidean norm in the (a;, v) phase space. Under suitable non-d egeneracy hypotheses 
for th e random flow defined by the dynamics, it was shown by iLedrappier fc Yound 
l)l988 t) that at small scales, rrir behaves as r''^ . More precisely, for almost every particle 
trajectory, 

ln^.(X(0,VW,0^^^ ^ 
Inr 

If the limit exists then it is equal to the information dimension of the statistical-steady- 
state density / (defined as the smallest Hausdorff dimension of a set of non-vanishing 
mass). In other words, the information dimension is equal to the Lyapunov dimension, 
identity which was conjectured by Kaplan & Yorke and which was rigorously proved by 
Ledrappier & Young for SRB measures. Of course, for a finite radius r of the ball, the 
mass of particles does not scale exactly as r'^^ but generally fluctuates. These fluctuations 
can be captured by investigating the scaling properties of the moments of order n of the 
mass, which follow power laws with exponents ^„ at small r's: 

(r<(X(t),V(t),<)> - r«" when r < i. (5.3) 

Here (•) denotes averaging with respect to the realisations of the carrier flow. These 
exponents can be expressed in term s of the spectrum of dimensions Z)„ fsee iGrassberged 
Il983t iHentschel fc Procaccialll983|) by = nDn+1- The result of Ledrappier fc Young 
clearly implies that (d^„ /dn) | n=o = d-L- Th ere is yet another exact result which relates the 
exponents ^„'s to the linearised dvnamics. iBaxendale fc Stroocld IfOSS.) showed that the 
exponent is given by t he Lyapunov mom ents dimension associated to the generalised 
Lyapunov exponents fsee lBenzi et a/.|ll985|) . This means that it satisfies 

lim i ln(|fi(i)r«0 =0, (5.4) 

t-^oo t 

where the vector R{t) denotes the phase-space separation between two infinitesimally 
close trajectorie s. Here we observe that phenomenological arguments are proposed in 
iBec et al\ ()2003(l to relate the full set of scaling exponents to the linearised system (|3.1|l 
governing the time evolution of the infinitesimal separation R{t). More precisely, for 
an arbitrary order n, the scaling exponent ^„ can be expressed in terms of the entropy 
function governing the long-time fluctuations of the stretching rates around their limit- 
ing values, the latter being precisely the Lyapunov exponents. We shall return to these 
matters and their implications in the conclusion (f©. 

Note that a self-similar distribution of particles would have ^„ — ucIl- Because of 
the convexity of the function n i— > ^„ (a consequence of a Holder's inequality applied 
to the moments of mass), it is clear that a necessary and sufficient condition for the 
exponents to depend nonlinearly on the order n is that < d^. Such a nonlinear 
behaviour is frequently referred to as anomalous scaling since it cannot be predicted by 
dimensional analysis. It implies that the particle mass distribution is intermittent and in 
particular that its probability density function (PDF) does not have Gaussian ta ils. More 
specifically, according to the usual multifractal formalism (see, e.g. . lFrisclilll995l) . we can 
write the stationary PDF of the phase-space mass m = mr(X(i), V(i), t) approximately 
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In?' 

The scaling (|5.3|l implies that in the limit r ^ with (lnm)/(lnr) — h fixed, the function 
D{r,h) tends to a limit that we denote DQ{h). The scaling exponents ^„ are related to 
Do{h) by a Legendre transform 

^„ = 2d + M[hn- Do{h)]. (5.6) 

h 

The limit DQ^h) is frequently referred to as the multifractal spectrum of the random 
dynamical attractor. The names entropy function, rate function or Cramer function as- 
sociated to the mass fluctuations are also sometimes used to refer to 2d — D^y{h). The 
multifractal spectrum can be interpreted heuristically as the dimension of the set on 
which the mass in small balls scales as r'* . A straightforward consequence of the almost- 
sure scaling of mass (|5.2(l is that the maximum of Do{h) is 2d and is attained for h = dh- 
Knowledge of Do{h) gives a very fine description of the small-scale properties of the mass 
distribution in phase space. 

5.2. Distribution of particle positions at low Stokes numbers 

In many applications and in particular those where interactions between inertial parti- 
cles are introduced (such as collisions or chemical reactions) , it is of greater interest to 
characterise not the phase-space mass but the position-space distribution of particles. To 
obtain the spatial density p{x, t) one must integrate the phase-space density / over the 
covelocities. As seen in ij4.2l when the dimension of the attractor is smaller than d, the 
density p is singular with support on the projection of the attractor, that is on clusters. 
To investigate the distribution of particles, we consider a small ball of radius r centred 
on a given position x which does not necessarily correspond to the position of a particle 
and hence may not intersect the attractor. The mass of particles contained in this ball is 



TO, 



{x,t)= dyp{x + y,t)= dy dv f{x,v,t). (5.7) 



Similarly to the phase-space mass distribution, the moments of to^ are expected to display 
anomalous scaling. The usual multifractal formalism used in the previous subsection leads 
to defining the multifractal spectrum Do(^) associated to the fluctuations of the position- 
space mass TOr- Relating the distribution of rur to that of the mass in phase space nir 
requires an integration over covelocities. This can generally not be done without a better 
understanding of the fractal properties of the attractor in the covelocity directions. It is 
clear that the position variables and the covelocity variables do not play similar roles in 
the dynamics, so that we have a kind of anisotropy in the phase-space dynamics. This 
prevents us from using standard tools of the geometry of fractal sets when integrating 
over the covelocities. 

This matter does however simplify if we limit ourselves to the asymptotics St <^ 1. 
The position-space mass distribution can then be obtained at intermediate spatial scales 
from the mass distribution in the full phase space. Indeed, the covelocities of the particles 
associated to low Stokes numbers are within a distance order St of the phase-space 
manifold v = {1 — j3) u{x,t). Hence, when St ^ r/L <C 1, the position-space mass rrir 
is given by the particles contained in a phase-space volume also of size r centred on 
this manifold (see figure . We thus introduce the exponents Cn for this intermediate 
asymptotic range: 

{Wir ") - for St < r/L < 1. (5.8) 
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{l-<f,)u(x,t) 




Figure 6. Sketch illustrating the arguments used at low Stokes number for the determination 
of the mass of particles rfir in a small volume of size r in the position space. The case d = 2 is 
here assumed and the two components of the covelocity in the four- dimensional phase space are 
schematically represented as a single variable. 



We can relate the Cn's to the exponents ^„ associated to the phase-space mass to^, taking 
into account that is evaluated along the attractor while rfir is not. We have seen in 21 
that for St -^l, the particle are located on dynamical clusters whose fractal dimension is 
equal to that of the phase-space attractor. The probability that the ball of size r intersects 
the fractal cluster of particles is hence cx r'^~'^H ^ ■\Yhen it intersects the projection of the 
attractor, the mass inside the ball obviously behaves as r^" . We thus have 

Cn ^d-d„+Cn- (5.9) 

The relation H5.2|) of Ledrappier & Young between the information dimension and the 
Lyapunov dimension implies that (dCn/dn)|„=o = dL- Moreover, the fact that = 0, 
together with relation 1)5. 9|l . implies that C,o = d — du- Conservation of the total mass of 
particles implies that the average mass in a ball of the position space scales as its volume. 
We hence have C,i = d and thus = dn- Note that because of the convexity of n i-^- 
a necessary and sufficient condition to have multifractal clusters of particles (that is a 
non-linear dependence of the C„) is that dn < dL and, of course, a value of the Stokes 
numb er below the c r itical clustering value. 

iHunt fc KaloshinI lll997l) showed that the dimension spectrum _D„ of fractal measures 
defined in the previous subsection is preserved under typical projections when 1 < n ^ 2 
("typical projections" means here almost all of them). From a direct application of this 
result to the steady-state phase-space density /, one expects that for order- unity values of 
the Stokes number, Cn = d— di/-|-min(nd, ^„) for < n < 1. When St is below the critical 
value, it is clear from convexity of the ^„'s together with (d^„/dn)|„=o = cLl < d that 
£,n < nd for all n. This implies that the relation 1)5. 9|l may be extended to finite values of 
St. However, it is not yet completely clear whether or not this result is applicable in our 
setting. Indeed, because of the aforementioned anisotropy of the phase-space dynamics, 
it is conceivable that the projection onto the position space is not a typical one. 

Our two-dimensional numerical experiments show that in the low Stokes number 
regime, the distribution of inertial particle positions has indeed multifractal properties. 
For St <^ 1, particle positions almost fill the whole space and, as seen previously, their 
Lyapunov dimension behaves as ~ d — C St^ . The scaling exponents Cn of the mo- 
ments of their mass distribution are hence very close to dn = 2 n, namely those given by 
simple tracers dynamics. Obtaining strong evidence for multifractal clustering requires 
an accurate determination of the exponents Cn; this in turn requires performing long av- 
erages over the realisations of the carrier flow. Figure shows the scaling exponents (^„'s 
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Figure 7. Scaling exponents C^„ of the moments of the mass of particles in a small ball of 
radius r as a function of their order n. The mass was calculated by box-counting with a total of 
A'^ = 10^ particles and for four different values of the Stokes number, as labelled. The solid line 
represents the exponents associated to a uniform distribution of particles. 
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Figure 8. Mass distribution for very heavy particles. The quantity D{r, h) defined by H5.10|l is 
plotted here for different non-dimensional radii r/L ranging from 5 X 10"^ to 2 x 10"^ and for 
two different values of the Stokes number as labelled. The density field is determined from the 
positions of A'' = 10^ particles. The curves approximately collapse on the function Do{h) which 
measures mass fluctuation at those scales. It attains its maximum dn for h — (11- 

obtained after time averaging the mass distribution for more than 10^ turnover times of 
the carrier flow. Here we take very heavy particles (/3 <C 1) and plot the scaling exponents 
C„ for different values of the Stokes number ranging from 10"'* to 10~*. Our results show 
that, already at very low Stokes number, the discrepancy from a space-filling distribution 
is increasingly noticeable when increasing the order n of the moment. 

The discrepancy from a uniform distribution is even clearer when measuring the fluc- 
tuations of the position-space mass distribution. In a way similar to what was done for 
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the multifractal analysis in phase space, we define 

D(^r,h)^d-^^^, (5.10) 
Inr 

where Pr{h) is the PDF of ft- = In m^/ Inr. Scaling of the position-space mass rfir in the 
intermediate asymptotics St -^r/L <^ 1 implies that in this range D{r, h) depends only 
weakly on the ball radius r and becomes thus essentially a function Do{h). It is easily 
checked that C„ = inih{n h + d — Do{h)). Together with H5.6|l and (|5.9(l this implies that 
Do{h) = du + Daih) — 2d. Thus, Do{h) has a maximum equal to dn for h — dL- 

In figure |S1 the functions D{r, h) corresponding to two different Stokes number are 
superposed for different values of the ball radius r in the intermediate asymptotics. 
These curves were obtained in the same setting as for figure|3 namely for infinitely heavy 
particles (/3 = 0) and by time averaging over 10^ turnover times the coarse-grained density 
obtained from the position of 10^ particles. Pretty good scaling is observed since the 
curves almost collapse on a same multifractal spectrum DQ{h). This implies in particular 
that the multifractal formalism catches well the mass distribution at scales within the 
intermediate range St <^ r / L <^ 1. 

The Lyapunov dimensions associated to the two small values of St investigated here 
are almost equal, from which one might infer that the mass distributions are the same. 
In fact this is not the case: the results of the numerical experiments presented here 
give two kinds of evidence for a different behaviour. First, for St <^ 1, the Hausdorff 
dimension of the attractor d}{ may not tend to d as fast as d^ (i.e. it may approach 
d slower than quadratically). This implies that clustering may be stronger than that 
predicted from the analysis of the Lyapunov exponents. Second, the large fluctuations 
of the mass distribution (both at large and low values) which are represented by the 
tails of the function Do{h) differ markedly for the two different low values of the Stokes 
numbers considered. At small Stokes numbers where the Stokes drag is very strong, the 
inertial particles almost follow the motion of the fluid. It is thus tempting to believe that 
the small effects due to particle inertia can be easily understood and quantified. As we 
have seen, this is not really the case: the small Stokes number asymptotic regime has a 
nontrivial clustering mechanism and there remains a number of open questions. 

6. Concluding remarks 

We have shown that the clustering of inertial particles in spatially smooth flow can 
be analysed using tools from dissipative dynamical systems. We thereby obtained evi- 
dence for the existence of a critical Stokes number below which particles form dynamical 
fractal clusters in position space, threshold whose existence for simple random flows was 
confirmed in both two and three dimensions by numerical simulations. We showed that 
the fluctuations in phase space of the singular mass distribution are derivable from the 
dimension spectrum of a random dynamical attractor. It is however difhcult to "project 
this down" to determining mass fluctuations merely in position space, because the inte- 
gration over velocities requires the knowledge of the fractal properties of the attractor 
in the phase-space directions associated to particle velocities. This cannot be done with- 
out a fuller understanding of the phase-space dynamics. So far, we have obtained only 
partial results in the asymptotics of low Stokes numbers, a regime with many practical 
applications. 

Let us mention some further questions arising in the low Stokes number limit. The 
evidence is that already at very low Stokes numbers, say 10"'^, particles form fractal 
clusters in position space. More precisely, our numerical results confirmed the prediction 
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of lBalkovskv et all ^OOl) that the discrepancy from a uniform distribution measured 
from the Lyapunov dimension goes quadraticaUy to zero as S'i — > 0. A delicate issue is 
whether or not the fractal Hausdorff dimension of the clusters tends to the dimension of 
the ambient space that fast. Obtaining the Hausdorff dimension numerically is difficult 
and requires extensive computational resources. Hence, a systematic investigation of its 
dependence on the dynamical parameters is at the moment out of reach. There is how- 
ever an alternative approach to some of the small Stokes number issues. As proposed by 
iMaxevI fl987). the particle dynamics at low Stokes numbers in an incompressible flow 
can be approximately described as an advection of simple passive tracers by a synthetic 
flow comprising a small compressible component (see also H4.2|l . In this approximation, 
we can then work directly in position space without having first to understand the full 
phase-space dynamics. In this simpler framework we can in principle take advantage of 
a recently proposed relation between the multifractal properties of the particle distri- 
bution (includin g the Hausdorff dimension) and the fluctuations of the stretching rates 
l|Bec a/.ll2003|) . For this we need to know the large-time statistical properties associ- 
ated to the local dynamics and, in particular, the large deviations of the stretching rates 
from their limiting values, the Lyapunov exponents. For the moment such a determina- 
tion has been made analytically only for spatially smooth compressible variants of the 
Kraichnan model, that is for velocity fields which are delta-correlated in time (see, e.g., 
iFalkovich et a l. 2001). The problem is that delta-correlated flow cannot be used for the 
low Stokes numbers regime. Indeed, the expansion requires that the Stokes numbers be 
much smaller than the (non-dimensional) correlation time of the carrier velocity field. 
A purely analytic investigation may thus be difficult, but we can use a mixed theoreti- 
cal and numerical strategy for determining the multifractal properties in the low Stokes 
number regime: use the aforementioned relation to express them in terms of the large 
deviations of the stretching rates, the latter being determined numerically. This should 
be much easier than direct numerical determination of multifractal exponents by, for 
example, box-counting methods. 

Finally, an important and non-trivial extension concerns the clustering properties of 
inertial particles at scales within the inertial range of turbulence. The dynamics of the 
flow at these scales is close to Kolmogorov's 1941 theory; thus the velocity is not smooth 
but, approximately, Holder continuous of exponent 1/3. As a consequence, tracer parti- 
cles separate in an explosive way, given by Richardson's t^^^ law. It is thus important 
to identify and to understand the mechanisms leading to preferential concentrations of 
particles at those scales where there is a strong competition between clustering due to 
dissipative dynamics and the explosive Richardson separation. This is a challenging ques- 
tions with serious methodological difficulties. On the one hand, numerical exploration of 
the properties of mass distribution at those scales may require very long time averages, 
as seen in the present study for smooth flows. Very long direct numerical simulations 
of three-dimensional high Reynolds number flows are hardly feasible. It may hence be 
necessary to focus on simpler flows, such as the synthetic turbulence of the Kraichnan 
model, the two-dimensional inverse cascade or large-eddy simulations. On the other hand, 
an analytical approach to this question can clearly not make use of standard dynamical 
systems tools which assume smooth dynamics. It thus seems of interest to adapt to the 
inertial particle dynamics models and tools which have been successfully applied in recent 
years to passive scalars advected by non-smooth flow in Kolmogorov-like regimes. 

I am deeply grateful to U. Frisch for many interesting discussions and for his continuous 
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